package picoevo.es.cmaes.cma;

import java.io.FileInputStream;
import java.io.FileNotFoundException;
import java.io.FileWriter;
import java.io.IOException;
import java.io.PrintStream;
import java.io.PrintWriter;
import java.util.Arrays;
import java.util.Comparator;
import java.util.Locale;
import java.util.Properties;
import java.util.Random;

import picoevo.es.cmaes.Jama.EigenvalueDecomposition;
import picoevo.es.cmaes.Jama.Matrix;

/* open questions:

 TODO: 

 Multistart is implemented in TestCMAExample1. This is more practical
 than with the CMAEvolutionStrategy class. 

 lambda and population size etc. needs to be cleaned up in class Parameters:
 setLambda() can be called before init(dim), all other will either not 
 have impact (setMu() goes through, but will be overwritten in init()). 
 setPopulationSize() cannot be called before init() (gives an error), 
 because this sets remaining parameters to their dim-dependent
 values. 

 Options: clean up? 

 // cma.setRecombination(cma.superlinear); TODO


 improve error management to reasonable standard 

 write and read object: 
 http://www.javacoffeebreak.com/articles/serialization/
 class ObjectOutputStream
 import java.io.*; 

 class SerializeTest implements Serializable{
 transient int x; 
 private int y; 
 public SerializeTest(int a, int b){
 x = a; 
 y = b; 
 }

 public String toString(){
 return "{x=" + x + ", y=" + y + "}"; 
 }

 public static void main(String[] args) throws Exception{

 SerializeTest st = new SerializeTest(66, 61); 
 System.out.println("Before Write := " + st);

 System.out.println("\n Writing SerializeTest object to disk");
 FileOutputStream out  = new FileOutputStream("serialized.txt");
 ObjectOutputStream so = new ObjectOutputStream(out);    
 so.writeObject(st);
 so.flush();

 System.out.println("\n Reading SerializeTest object from disk\n");
 FileInputStream in     = new FileInputStream("serialized.txt");
 ObjectInputStream si   = new ObjectInputStream(in); 
 SerializeTest fromdisk = (SerializeTest)si.readObject();

 /* x will be 0 because it won't be read from disk since transient 
 System.out.println("After Read := " + fromdisk);
 }
 }

 See TODO in source code

 DONE: 

 StopCondition: should it really be an inner class? Yes, because it
 heavily uses the class CMAEvolutionStrategy

 */

//import Jama.*;
//import Jampack.*;
//import weka.core.*;
/**
 * implements an Evolution Strategy with Covariance Matrix Adaptation for
 * non-linear, non-convex global function minization. The CMA-Evolution Strategy
 * (CMA-ES) is a reliable stochastic optimization method which should be
 * applied, if derivative based methods, e.g. quasi-Newton BFGS or conjucate
 * gradient, fail due to a rugged search landscape (e.g. noise, local optima,
 * outlier, etc.) of the objective function. Like a quasi-Newton method the
 * CMA-ES learns and applies a variable metric of the underlying search space.
 * Unlike a quasi-Newton method the CMA-ES does neither estimate nor use
 * gradients making it considerably more reliable (in terms of finally finding a
 * good, or close to optimal, solution).
 * 
 * <p>
 * In general, on smooth objective functions the CMA-ES is roughly ten times
 * slower than BFGS (counting objective function evaluations, no gradients
 * provided). For up to dim=10 variables even the derivative-free simplex direct
 * search method (Nelder & Mead) can be faster, but it is far less reliable than
 * CMA-ES.
 * 
 * <p>
 * The CMA-ES is particularely well suited for non-separable and/or badly
 * conditioned problems. To observe the advantage of CMA, compared to a
 * conventional evolution strategy, it will usually take about 30&#215;dim
 * function evaluations. On difficult problems the complete optimization (a
 * single run) is expected to take <em>roughly</em> at least
 * 20&#215;dim&#215;min(10,dim) and up to 300&#215;(dim+2)<sup>2</sup>
 * function evaluations.
 * 
 * <p>
 * There are four ways to use the CMAEvolutionStrategy class.
 * <ol>
 * <li> via the samplePopulation() and updateDistribution() methods
 * (preferably).
 * <li> via the sampleSingle() and updateSingle() methods.
 * <li> via the iterate() method using the IObjectiveFunctionParallel interface.
 * <li> via the evolveSingle() method using the IObjectiveFunction interface.
 * </ol>
 * In general, these ways cannot be mixed. If the objective function can be
 * evaluated in parallel the methods samplePopulation() or iterate() should be
 * used.
 * 
 * <p>
 * Example code snippet:
 * 
 * <pre>
 *   CMAEvolutionStrategy cma = new CMAEvolutionStrategy();
 *   cma.loadProperties(&quot;filename&quot;); // read options, see file CMAEvolutionStrategy.properties
 *   cma.setDimension(10);           // overwrite some loaded properties
 *   cma.setInitialSearchRegion(-1, 1);
 *   cma.opts.stopFunctionValue = new Double(1e-6); 
 *   while(cma.stopConditions.isFalse()) {
 *       double[][] pop = cma.samplePopulation();
 *       double[] functionValues;
 *       for (int i = 0; i &lt; pop.length; ++i) {
 *          ...assign 
 *          functionValues[i] = objective function value of pop[i] here...
 *       }
 *       cma.updateDistribution(functionValues);
 *       cma.writeToDefaultFilesNew();
 *       ...in case, print output, eg. System.out.println(cma.print()), 
 *          or process best found solution, getBestSolution()...
 *   } // while 
 *   ...retrieve best solution, termination criterion via stopConditions etc...    
 * </pre>
 * 
 * Another code example class is TestCMAExample1.
 * 
 * 
 * @author Nikolaus Hansen, 1996, 2003, 2005
 */
public class CMAEvolutionStrategy implements java.io.Serializable {
	/**
	 * 
	 */
	private static final long serialVersionUID = -7237782648003486603L;
	public final String versionNumber = new String("0.91");

	// [!n] stable mode. Introduced by Nicolas B. to override restrictions on
	// setting boundaries (Niko Hansen says that current implementation is not
	// optimal - override this flag at your own risks)
	private boolean _stableMode = true;

	/**
	 * Interface to (mostly generic) options for the optimization, like the
	 * maximum number of objective function evaluations. Options are
	 * <em>not</em> needed to initialize CMAEvolutionStrategy and can be set
	 * and changed any time. For the moment the interface works by directly
	 * assigning new values to the fields.
	 */
	public class Options { // has full access to the embracing class
		// needs to be public to make sure that a using class can excess
		// Options.
		// Therefore, if not nested, needs to move into a seperate file
		// TODO: rethink how to use these options inside and outside,
		// it should become a "real class" with setter methods for all fields!?

		/**
		 * lower bound for standard deviations (step sizes). The Array can be of
		 * any length. The i-th entry corresponds to the i-th variable. If
		 * length&#60;dim the last entry is recycled for all remaining
		 * variables. Zero entries mean, naturally, no lower bound. Be aware of
		 * the interference with stopTolX (and stopTolXFactor): if lowerStdDev
		 * is larger than stopTolX, the termination criterion can never be
		 * satisfied.
		 * <p>
		 * Example:
		 * 
		 * <pre>
		 * lowerStdDev = new double[] { 1e-4, 1e-8 };
		 * </pre>
		 */
		public double[] lowerStandardDeviations;
		/**
		 * upper bound for standard deviations (step lengths). Zero entries mean
		 * no upper bound. Be aware of the interference with stopTolUpXFactor.
		 * 
		 * @see #lowerStandardDeviations
		 * @see #stopTolUpXFactor
		 */
		public double[] upperStandardDeviations;
		/**
		 * stop if function value drops below the target function value
		 * stopFunctionValue. Default = -inf.
		 */
		public Double stopFunctionValue; // only used if non-null
		/**
		 * stop if function value differences within a single generation become
		 * smaller than stopTolFun. Default = 1e-12. The measured "function
		 * value differences" include repairment penalties (eg. for boundary
		 * repairment).
		 */
		public Double stopTolFun = new Double(1e-12); // only used if non-null
		/**
		 * stop if function value differences within the last 10 + 30*lambda/N
		 * generations become smaller than stopTolFunHist. Default = 1e-12. The
		 * measured objective function value differences do not include
		 * repairment penalties.
		 */
		public Double stopTolFunHist = new Double(1e-12); // used if non-null
		/** stop if search steps become smaller than stopTolX. Default = 0 */
		public Double stopTolX = new Double(0.0); // used if non-null
		/**
		 * stop if search steps become smaller than stopTolXFactor * initial
		 * step size. Default = 1e-9.
		 */
		public Double stopTolXFactor = new Double(1e-9); // used if TolX is
		// null
		/**
		 * stop if search steps become larger than stopTolUpXFactor * initial
		 * step size. Default = 1e3. When this termination criterion applies the
		 * initial step-size is chosen far too small or divergent behavior of
		 * the CMA-ES is observed.
		 */
		public Double stopTolUpXFactor = new Double(1e3); // multiplier for
		// initial sigma
		/**
		 * stop if the number of objective function evaluations exceed
		 * stopEvaluations
		 */
		public Long stopEvaluations = new Long(Long.MAX_VALUE); // TODO: set
		// default when
		// dimension is
		// defined
		/** stop if the number of iterations (generations) exceed stopIterations */
		public Long stopIterations = new Long(Long.MAX_VALUE);
		/**
		 * flag used by methods iterate() and evolveSingle(), whether to write
		 * output to files. Methods write an output file if flgWriteFile&#62;0.
		 * 
		 * @see #write()
		 */
		public int flgWriteFile = 0;
		/**
		 * determines whether writeTo...-methods overwrite output files
		 * initially
		 */
		public boolean flgAppendFiles = false;
		/** Output files written will have the names outputFileNamesPrefix*.dat */
		public String outputFileNamesPrefix = "cmaes";
		/** verbosity in terms of printing and writing files, not yet clear cut */
		public int verbose = 2;
		public Double maximumTimeFractionForEigendecomposition = new Double(0.5);

		void setOptions(Properties properties) {
			String s;
			if ((s = properties.getProperty("stopFunctionValue")) != null)
				stopFunctionValue = Double.valueOf(getFirstToken(s));
			stopTolFun = getFirstToken(properties.getProperty("stopTolFun"), stopTolFun);
			stopTolFunHist = getFirstToken(properties.getProperty("stopTolFunHist"), stopTolFunHist);
			stopTolX = getFirstToken(properties.getProperty("stopTolX"), stopTolX);
			stopTolXFactor = getFirstToken(properties.getProperty("stopTolXFactor"), stopTolXFactor);
			stopTolUpXFactor = getFirstToken(properties.getProperty("stopTolUpXFactor"), stopTolUpXFactor);
			stopEvaluations = getFirstToken(properties.getProperty("stopEvaluations"), stopEvaluations);
			stopIterations = getFirstToken(properties.getProperty("stopIterations"), stopIterations);
			if ((s = properties.getProperty("upperStandardDeviations")) != null && !s.equals(""))
				upperStandardDeviations = parseDouble(getAllToken(s));
			if ((s = properties.getProperty("lowerStandardDeviations")) != null && !s.equals(""))
				lowerStandardDeviations = parseDouble(getAllToken(s));
			outputFileNamesPrefix = properties.getProperty("outputFileNamesPrefix", outputFileNamesPrefix).split("\\s")[0];
			maximumTimeFractionForEigendecomposition = getFirstToken(properties.getProperty("maximumTimeFractionForEigendecomposition"), maximumTimeFractionForEigendecomposition);
			// stop = Double.valueOf(properties.getProperty("stop",
			// Double.toString(stop.doubleValue())));

		}
	}

	/**
	 * Interface to strategy parameters for the CMA Evolution Strategy. The
	 * class CMAEvolutionStrategy.Parameters processes the strategy parameters,
	 * like population size and learning rates, for the class
	 * CMAEvolutionStrategy. The main issues are dependencies (eg, the parent
	 * number, mu, cannot be larger than the population size lambda) and
	 * respective consistency checking. Set methods are available, but not all
	 * parameters can be changed after the first solution search point is
	 * delivered via CMA methods sample*().
	 * 
	 * <p>
	 * Most commonly, the offspring population size lambda can be changed
	 * (increased) from its default value to improve the global search
	 * capability via setPopulationSize. We recommend to first use the default
	 * values.
	 * 
	 * @see #setPopulationSize(int)
	 */
	public class Parameters {
		// int N; // uses N from outer class, see also init()
		int lambda;
		int mu;
		double[] weights;
		double mueff;
		double mucov;
		double ccov = -1;
		double damps;
		double cc;
		double cs;

		double chiN;

		int status = -1;

		public Parameters() {
		}

		public Parameters(int N) {
			init(N);
		}

		/**
		 * Initialization of parameters with defaults depending on lambda
		 * (offspring population size) that might have been set already.
		 * 
		 * @param dimension
		 *            integer search space dimension, number of variables
		 */
		public void init(int dimension) {
			// this.N = N;
			N = dimension;
			chiN = Math.sqrt(N) * (1.0 - 1.0 / (4.0 * N) + 1.0 / (21.0 * N * N));
			setDefaults();
		}

		public void setDefaults() {
			if (N == 0)
				error("cannot set defaults without knowing dimension");
			if (lambda == 0)
				setLambda((int) (4.0 + 3.0 * Math.log(N)));
			setDefaultsFromPopsize();
		}

		public void setDefaultsFromPopsize() {
			setMu();
			setWeights(getMu());
			setMucov();
			setCcov();
			setDamps();
			setCs();
			// setCc();
			makeValid();
			status = 2;
		}

		public void makeValid() {
			if (getMu() > getLambda()) {
				setLambda(2 * getMu());
			}
			if (getMu() != getWeights().length) {
				setWeights(getMu());
			}
			if (getLambda() < 5)
				setLambda(5);
			setCcov(ccov);
			setCc(cc);
			setCs(cs);
			setDamps(damps);
			status = 1;
		}

		/**
		 * Getter for property mu.
		 * 
		 * @return Value of property mu.
		 * 
		 */
		public int getMu() {
			if (mu == 0)
				setMu();
			return mu;
		}

		/**
		 * Sets default value for parent number mu.
		 * 
		 * @see #setMu(int mu)
		 * @see #setSelectionFraction
		 */
		public void setMu() {
			setMu(getLambda() / 2);
		}

		/**
		 * Setter for parent number mu, use preferably setSelectionFraction.
		 * 
		 * @param mu
		 *            New value for the number of parents mu.
		 * @see #setSelectionFraction
		 */
		public void setMu(int mu) {
			this.mu = mu;
			if (weights == null)
				setWeights(mu);
			else if (weights.length != mu)
				setWeights(mu);
		}

		/**
		 * Getter for offspring population size lambda
		 * 
		 * @return Value of lambda
		 * 
		 */
		public int getLambda() {
			if (lambda == 0)
				setLambda();
			return lambda;
		}

		/**
		 * Setter for population size alias sample size to default value
		 * (int)(4.0+3.0*Math.log(N))
		 */
		protected void setLambda() {
			if (N == 0)
				error("cannot set lambda to default without knowing dimension");
			setLambda((int) (4.0 + 3.0 * Math.log(N)));
		}

		int flgLambdaChanged = 0;

		/**
		 * Setter for offspring population size alias sample size alias lambda
		 * 
		 * @param lambda
		 *            set population size
		 */
		protected void setLambda(int lambda) {
			if (arx != null)
				if (arx.length != 0 && arx.length != lambda) {
					flgLambdaChanged = 1;
					if (1 < 3)
						error("population size cannot be changed at this point");
				}
			// TODO: make population size dynamic
			// effected are arx, arxrepaired,...??...
			// only at state==0 and 3 sensible?
			this.lambda = lambda;
		}

		/** @see #getLambda() */
		public int getPopulationSize() {
			return getLambda();
		}

		/**
		 * Setter for offspring population size (lambda) to default value
		 * (int)(4.0+3.0*Math.log(N)) and depending parameters, like parent
		 * number (mu), accordingly
		 */
		public void setPopulationSize() {
			setLambda();
			setDefaultsFromPopsize();
		}

		/**
		 * Setter for offspring population size (lambda) and parent number (mu)
		 * accordingly
		 * 
		 * @param lambda
		 *            is the offspring population size
		 */
		public void setPopulationSize(int lambda) {
			setLambda(lambda);
			setDefaultsFromPopsize();
		}

		public final int superlinearRecombinationWeights = 0;
		public final int linearRecombinationWeights = 1;
		public final int equalRecombinationWeights = 2;

		/**
		 * TODO: to be implemented. Set the fraction of the offspring population
		 * sample selected for updating the distribution for the next generation
		 * and the method to choose the recombination. The parent number mu is
		 * approximately fraction&#215;lambda, where lambda is the offspring
		 * population size.
		 * 
		 * @param fraction
		 */
		public void setSelectionFraction(double fraction, int recombinationType) {

		}

		/**
		 * Getter for property weights.
		 * 
		 * @return Value of property weights.
		 * 
		 */
		public double[] getWeights() {
			if (weights.length == 0)
				setWeights(getMu());
			return this.weights;
		}

		/**
		 * TODO: to be implemented. Recombination weights can be equal, linear
		 * decreasing, or superlinear decreasing. The respective parameter is
		 * CMAEvolutionStrategy.Parameters.equalRecombinationWeights,
		 * ...linear..., or ...superlinear...(default) The parent number mu will
		 * be set accordingly to 0.25 lambda for equal weights, and 0.5 lambda
		 * otherwise.
		 * 
		 * @param recombinationType
		 * @see #setSelectionFraction
		 */
		public void setRecombinationWeights(int recombinationType) {
			setWeights(getMu());
		}

		/**
		 * Setter for recombination weights
		 * 
		 * @param mu
		 *            is the number of parents, number of weights > 0
		 */
		private void setWeights(int mu) {
			double[] weights = new double[mu];
			for (int i = 0; i < mu; ++i) {
				weights[i] = (Math.log(mu + 1) - Math.log(i + 1));
			}
			setWeights(weights);
		}

		protected void setWeights(double[] weights) {
			double sum = 0;
			for (int i = 0; i < weights.length; ++i)
				sum += weights[i];
			for (int i = 0; i < weights.length; ++i)
				weights[i] /= sum;
			this.weights = weights;
			setMu(weights.length);
			setMueff(weights);
		}

		/**
		 * Getter for property mueff, the "variance effective selection mass".
		 * 
		 * @return Value of property mueff.
		 * 
		 */
		public double getMueff() {
			if (mueff == 0)
				setMueff();
			return mueff;
		}

		/**
		 * Setter for property mueff.
		 * 
		 * @param mueff
		 *            New value of property mueff.
		 * 
		 */
		private void setMueff() {
			setMueff(getWeights());
		}

		private void setMueff(double[] weights) {
			double sum1 = 0;
			double sum2 = 0;
			for (int i = 0; i < mu; ++i) {
				sum1 += weights[i];
				sum2 += weights[i] * weights[i];
			}
			this.mueff = sum1 * sum1 / sum2;
			setMucov(this.mueff);
		}

		/**
		 * Getter for property mucov. mucov determines the mixing between
		 * rank-one and rank-mu update. For mucov = 1, no rank-mu updated takes
		 * place.
		 * 
		 * @return Value of property mucov.
		 * 
		 */
		public double getMucov() {
			if (mucov == 0)
				setMucov();
			return mucov;
		}

		/**
		 * Setter for mucov to the default mueff.
		 * 
		 * @see #getMucov()
		 */
		public void setMucov() {
			setMucov(getMueff());
		}

		/**
		 * Setter for mucov.
		 * 
		 * @param mucov
		 *            New value of mucov.
		 * @see #getMucov()
		 */
		public void setMucov(double mucov) {
			this.mucov = mucov;
		}

		/**
		 * Getter for property covariance matrix learning rate ccov.
		 * 
		 * @return Value of property ccov.
		 * 
		 */
		public double getCcov() {
			if (N == 0)
				warning("return value for ccov is questionable");
			return ccov;
		}

		/**
		 * Setter for covariance matrix learning rate ccov to default value.
		 * 
		 * @see #getCcov()
		 */
		public void setCcov() {
			if (N == 0)
				error("cannot set ccov to default without knowing dimension");
			// double ccov = 2.0 / (N + 1.41) / (N + 1.41);
			setCcov(2.0 / (N + 1.41) / (N + 1.41) / getMucov() + (1 - (1.0 / getMucov())) * Math.min(1, (2 * getMueff() - 1) / (getMueff() + (N + 2) * (N + 2))));
		}

		/**
		 * Setter for covariance matrix learning rate ccov.
		 * 
		 * @param ccov
		 *            New value of property ccov.
		 * @see #getCcov()
		 */
		public void setCcov(double ccov) {
			if (ccov < 0 || ccov > 1) {
				warning("setCcov(ccov) failed due to invalid value ccov=" + new Double(ccov).toString());
				setCcov();
			}
			this.ccov = ccov;
		}

		/**
		 * Getter for step-size damping damps. The damping damps determines the
		 * amount of step size change.
		 * 
		 * @return Value of damps.
		 * 
		 */
		public double getDamps() {
			if (damps <= 0)
				setDamps();
			return damps;
		}

		/**
		 * Setter for damps to default value.
		 */
		public void setDamps() {
			if (N == 0)
				error("cannot set damps to default without knowing dimension");

			// setDamps(Math.max(1.0, 3.0*getMueff()/(N+10.)) + getCs());
			setDamps(1 + 2 * Math.max(0, Math.sqrt((getMueff() - 1.) / (N + 1.)) - 1) + getCs());
		}

		/**
		 * Setter for damps.
		 * 
		 * @param damps
		 *            New value of damps.
		 * @see #getDamps()
		 */
		public void setDamps(double damps) {
			this.damps = damps;
		}

		/**
		 * Getter for backward time horizon parameter cc for distribution
		 * cumulation (for evolution path p<sub>c</sub>).
		 * 
		 * @return Value of cc.
		 * 
		 */
		public double getCc() {
			if (cc == 0)
				setCc();
			return cc;
		}

		/**
		 * Setter for cc to default value.
		 * 
		 */
		public void setCc() {
			if (N == 0)
				error("cannot set cc to default without knowing dimension");
			setCc(4.0 / (N + 4.0));
		}

		public void setCc(double cc) {
			if (cc <= 0 || cc > 1)
				setCc();
			this.cc = cc;
		}

		/**
		 * Getter for cs, parameter for the backward time horizon for the
		 * cumulation for sigma.
		 * 
		 * @return Value of property cs.
		 * 
		 */
		public double getCs() {
			if (cs == 0)
				setCs();
			return cs;
		}

		/**
		 * Setter for cs to default value.
		 * 
		 * @see #getCs()
		 */
		public void setCs() {
			if (N == 0)
				error("cannot set cc to default without knowing dimension");
			setCs((getMueff() + 2.0) / (N + (getMueff() + 3.0)));
		}

		/**
		 * Setter for cs.
		 * 
		 * @see #getCs()
		 */
		public void setCs(double cs) {
			if (cs <= 0 || cs > 1)
				setCs();
			this.cs = cs;
		}

	}

	/**
	 * Interface to termination criteria and status.
	 */
	public class StopCondition {
		int index = 0; // number of messages collected == index where to write
		// next message
		String[] messages; // Initialisierung
		double lastcounteval;

		/**
		 * true whenever a termination criterion was met. clear() re-sets this
		 * value to false.
		 * 
		 * @see #clear()
		 */
		public boolean isTrue() {
			return test() > 0 ? true : false;
		}

		/**
		 * evaluates to NOT isTrue().
		 * 
		 * @see #isTrue()
		 */
		public boolean isFalse() {
			return !isTrue();
		}

		/**
		 * greater than zero whenever a termination criterion was satisfied.
		 * clear() re-sets this value to zero.
		 * 
		 * @return number of generated termination condition messages
		 */
		public int getNumber() {
			return test();
		}

		/**
		 * get description messages of satisfied termination criteria. The
		 * messages start with one of "FunctionValue:", "TolFun:",
		 * "TolFunHist:", "TolX:", "TolUpX:", "Evaluations:", "Iterations:",
		 * "ConditionNumber:", "NoEffectAxis:", "NoEffectCoordinate:".
		 * 
		 * @return String[] s with messages of termination conditions.
		 *         s[0].equals("") is true if no termination condition is
		 *         satisfied yet
		 */
		public String[] getStrings() {
			test();
			if (index > 0)
				return messages; /* first string might be empty */
			else
				return new String[] { "" }; // TODO: decide about interface
			// before release
		}

		/**
		 * remove all earlier termination condition messages
		 */
		public void clear() {
			messages = new String[] { "" };
			index = 0;
		}

		private void append(String s) {
			// TODO could be replaced by
			// type ArrayList<String> or Vector<String>
			// disadvantage: parameterized types are not available before
			// Version 5.0
			String[] oldm = messages;
			messages = new String[index + 1];
			messages[index] = s + " (iter=" + new Long(countiter).toString() + ",eval=" + new Long(counteval).toString() + ")";
			/* copy old messages */
			for (int i = 0; i < index; ++i)
				messages[i] = oldm[i];
			++index;
		}

		/**
		 * Tests termination criteria and evaluates to non-zero when a
		 * termination criterion is satisfied
		 * 
		 * @return non-zero int whenever a termination criterion is satisfied
		 */
		int test() {
			if (state < 0)
				return 0; // not yet initialized
			if (state > 0 && lastcounteval == counteval)
				return index;
			lastcounteval = counteval;

			/* FUNCTION VALUE */
			if (opts.stopFunctionValue != null && (countiter > 1 || (countiter == 1 && state > 2))) {
				if (bestever.getFunctionValue() <= opts.stopFunctionValue.doubleValue())
					append("FunctionValue: Objective function value dropped below the target function value " + opts.stopFunctionValue);
			}

			/* #Fevals */
			if (counteval >= opts.stopEvaluations.longValue())
				append("Evaluations: maximum number of function evaluations " + opts.stopEvaluations + " reached");

			/* #iterations */
			if (countiter >= opts.stopIterations.longValue())
				append("Iterations: maximum number of iterations reached");

			/* TOLFUN */
			if (opts.stopTolFun != null && (countiter > 1 || (countiter == 1 && state >= 3))) {
				double tolfun = opts.stopTolFun.doubleValue();
				if (fit.fitness[fit.fitness.length - 1].d - fit.fitness[0].d <= tolfun)
					append("TolFun: function value changes below stopTolFun=" + tolfun);
			}

			/* TOLFUNHIST */
			if (opts.stopTolFunHist != null && countiter > fit.history.length) {
				double tolfun = opts.stopTolFunHist.doubleValue();
				if (math.max(fit.history) - math.min(fit.history) <= tolfun)
					append("TolFunHist: history of function value changes below stopTolFunHist=" + tolfun);
			}

			/* TOLX */
			double tolx = opts.stopTolXFactor.doubleValue() * minstartsigma;
			if (opts.stopTolX != null)
				tolx = Math.max(tolx, opts.stopTolX.doubleValue());
			if (sigma * maxsqrtdiagC < tolx && sigma * math.max(math.abs(pc)) < tolx)
				append("TolX: standard deviation below stopTolX=" + tolx);

			/* TOLXUP */
			if (sigma * maxsqrtdiagC > opts.stopTolUpXFactor.doubleValue() * maxstartsigma)
				append("TolUpX: standard deviation increased by stopTolUpXFactor=" + opts.stopTolUpXFactor.doubleValue() * maxstartsigma);

			/* Internal (numerical) stopping termination criteria */

			/* Jede Hauptachse i testen, ob x == x + 0.1 * sigma * rgD[i] * B[i] */
			for (int iAchse = 0; iAchse < N; ++iAchse) {
				int iKoo;
				double fac = 0.1 * sigma * diagD[iAchse];
				for (iKoo = 0; iKoo < N; ++iKoo) {
					if (xmean[iKoo] != xmean[iKoo] + fac * B[iKoo][iAchse])
						break; // is
					// OK
					// for
					// this
					// iAchse
				}
				if (iKoo == N) // no break, therefore no change for axis iAchse
					append("NoEffectAxis: Mutation " + new Double(0.1 * sigma * diagD[iAchse]).toString() + " in a principal axis direction has no effect");
			} /* for iAchse */

			/* Testen ob eine Komponente des Objektparameters festhaengt */
			for (int iKoo = 0; iKoo < N; ++iKoo) {
				if (xmean[iKoo] == xmean[iKoo] + 0.2 * sigma * Math.sqrt(C[iKoo][iKoo]))
					append("NoEffectCoordinate: Mutation of size " + new Double(0.2 * sigma * Math.sqrt(C[iKoo][iKoo])).toString() + "\n  in coordinate " + new Integer(iKoo).toString() + " has no effect");
			} /* for iKoo */

			/* Condition number */
			if (math.min(diagD) == 0)
				append("ConditionNumber: smallest eigenvalue equals zero");
			else if (math.max(diagD) / math.min(diagD) > 1e7)
				append("ConditionNumber: condition number of the covariance matrix exceeds 1e14");

			return index;
		}
	} // StopCondtion

	void testAndCorrectNumerics() {
		/* Test if function values are identical */
		if (getCountIter() > 1 || (getCountIter() == 1 && state >= 3))
			if (fit.fitness[0].d == fit.fitness[1 + sp.getLambda() / 2].d) {
				warning("flat fitness landscape, step-size increased");
				sigma *= Math.exp(0.2 + sp.getCs() / sp.getDamps());
			}
		// if (math.max(fit.history) == math.min(fit.history))
		// sigma *= Math.exp(0.2+sp.getCs()/sp.getDamps()); // or better
		// stop<->option?

		/* Align (renormalize) scale C (and consequently sigma) */
		/*
		 * e.g. for infinite stationary state simulations (noise handling needs
		 * to be introduced for that)
		 */
		double fac = 1;
		if (math.min(diagD) < 1e-9)
			fac = 1. / math.max(diagD);
		else if (math.max(diagD) > 1e4)
			fac = 1. / math.max(diagD);

		if (fac != 1.) {
			sigma /= fac;
			for (int i = 0; i < N; ++i) {
				pc[i] *= fac;
				diagD[i] *= fac;
				for (int j = 0; j < N; ++j)
					C[i][j] *= fac * fac;
			}
		}
	} // Test...

	/**
	 * fields of opts can be changed at any time to control the running behavior
	 */
	public Options opts = new Options();

	private Parameters sp = new Parameters(); // alias for inside use
	/**
	 * strategy parameters that can be set initially
	 */
	public Parameters parameters = sp; // for outside use also

	/** permits access to whether and which termination conditions were satisfied */
	public StopCondition stopConditions = new StopCondition();

	int N;
	long seed = System.currentTimeMillis();
	Random rand = new Random(seed); // Note: it also Serializable

	final MyMath math = new MyMath();
	double axisratio;
	public long counteval; // [n]
	long countiter;

	SolutionGeneric bestever;

	double sigma = 0.0;
	double[] initXLow, initXUp; // initial search
	double[] LBound, UBound; // bounds
	double[] xmean;
	double[] pc;
	double[] ps;
	Matrix mC;
	double[][] C;
	double maxsqrtdiagC;
	double minsqrtdiagC;
	Matrix mB;
	double[][] B;
	double[] diagD;

	/* init information */
	double[] startsigma;
	double maxstartsigma;
	double minstartsigma;

	boolean iniphase;

	/**
	 * state (postconditions): -1 not yet initialized 0 initialized init() 0.5
	 * reSizePopulation 1 samplePopulation, sampleSingle, reSampleSingle 2.5
	 * updateSingle 3 updateDistribution
	 */
	double state = -1;
	int mode = 0;
	final int SINGLE_MODE = 1;
	final int PARALLEL_MODE = 2;

	long countCupdatesSinceEigenupdate;

	/* fitness information */
	class FitnessCollector {
		double history[];
		IntDouble[] fitness; // int holds index for respective arx
		IntDouble[] funValues; // sorted differently than fitness!
		/**
		 * history of delta fitness / sigma^2. Here delta fitness is the minimum
		 * of fitness value differences with distance lambda/2 in the ranking.
		 */
		double[] deltaFitHist = new double[5];
		int idxDeltaFitHist = 0;
	}

	protected FitnessCollector fit = new FitnessCollector();

	double recentFunctionValue;
	double recentMaxFunctionValue;
	double recentMinFunctionValue;
	int idxRecentOffspring;

	double[][] arx;
	double[][] arxrepaired;
	double[] xold;

	double[] BDz;
	double[] artmp;

	// OptiPoint[] arpoints;

	String propertiesFileName = new String("CMAEvolutionStrategy.properties");

	/**
	 * postpones most initialization. For initialization use setInitial...
	 * methods or set up a properties file, see file
	 * "CMAEvolutionStrategy.properties".
	 */
	public CMAEvolutionStrategy() {
		state = -1;
	}

	/**
	 * initialization from properties input, see file
	 * "CMAEvolutionStrategy.properties" for valid properties
	 */
	public CMAEvolutionStrategy(Properties properties) {
		setProperties(properties);
		state = -1;
	}

	/**
	 * reads properties (options, strategy parameter settings) from file, see
	 * file "CMAEvolutionStrategy.properties", use "" for no property file
	 */
	public CMAEvolutionStrategy(String propertiesFileName) {
		this.propertiesFileName = propertiesFileName;
		state = -1;
	}

	/**
	 * @param dimension
	 *            search space dimension, dimension of the objective functions
	 *            preimage, number of variables
	 */
	public CMAEvolutionStrategy(int dimension) {
		setDimension(dimension);
		state = -1;
	}

	/**
	 * Forces setup of everything, is called once in samplePopulation(), that
	 * is, as late as possible. Can be, but needs not be to be called from
	 * outside
	 */
	public void init() {
		int i;
		if (N <= 0)
			error("dimension needs to be determined, use eg. setDimension() or setInitialX()");

		sp.init(N);
		sp = parameters; /* just in case the user assigned parameters */
		diagD = new double[N];
		for (i = 0; i < N; ++i)
			diagD[i] = 1;

		/* expand initial search region */
		if (initXLow != null && initXUp != null) {
			if (initXLow.length == 1) {
				double d = initXLow[0];
				initXLow = new double[N];
				for (i = 0; i < N; ++i)
					initXLow[i] = d;
			}
			if (initXUp.length == 1) {
				double d = initXUp[0];
				initXUp = new double[N];
				for (i = 0; i < N; ++i)
					initXUp[i] = d;
			}
		}

		/* expand Boundaries */
		if (LBound == null) {
			LBound = new double[N];
			for (i = 0; i < N; ++i)
				LBound[i] = Double.NEGATIVE_INFINITY;
		} else if (LBound.length == 1) {
			double d = LBound[0];
			LBound = new double[N];
			for (i = 0; i < N; ++i)
				LBound[i] = d;
		}
		if (UBound == null) {
			UBound = new double[N];
			for (i = 0; i < N; ++i)
				UBound[i] = Double.POSITIVE_INFINITY;
		} else if (UBound.length == 1) {
			double d = UBound[0];
			UBound = new double[N];
			for (i = 0; i < N; ++i)
				UBound[i] = d;
		}

		/* Initialization of xmean */
		/* expand if already set */
		if (xmean != null) {
			if (xmean.length == 1) {
				double d = xmean[0];
				xmean = new double[N];
				for (i = 0; i < N; ++i)
					xmean[i] = d;
			}
			/* set via initial search region */
		} else if (initXLow != null && initXUp != null) {
			xmean = new double[N];
			for (i = 0; i < N; ++i)
				xmean[i] = initXLow[i] + rand.nextDouble() * (initXUp[i] - initXLow[i]);
			/* set via boundaries */
		} else if (math.max(UBound) < Double.MAX_VALUE && math.min(LBound) > -Double.MAX_VALUE) {
			xmean = new double[N];
			for (i = 0; i < N; ++i)
				xmean[i] = LBound[i] + rand.nextDouble() * (UBound[i] - LBound[i]);
		} else {
			error("no initial search point (solution) X or search interval specified");
			xmean = new double[N];
			for (i = 0; i < N; ++i)
				xmean[i] = rand.nextDouble();
		}

		/* Initialization of sigmas */
		if (startsigma != null) { // 
			if (startsigma.length == 1) {
				sigma = startsigma[0];
			} else if (startsigma.length == N) {
				sigma = math.max(startsigma);
				if (sigma <= 0)
					error("initial standard deviation sigma must be positive");
				for (i = 0; i < N; ++i) {
					diagD[i] = startsigma[i] / sigma;
				}
			} else
				error("BUG detected");
		} else if (initXLow != null && initXUp != null) { // using initial
															// search
			// region
			for (i = 0; i < N; ++i)
				diagD[i] = 0.3 * (initXUp[i] - initXLow[i]);
			if (math.min(diagD) <= 0) { /* don't set sigma to zero */
				error("initial search region has zero volume, use eg. setInitialStandardDeviations() to set 'sigma'");
				for (i = 0; i < N; ++i)
					diagD[i] = 1.;
			} else {
				sigma = math.max(diagD);
				if (sigma <= 0)
					error("initial standard deviation sigma must be positive");
				for (i = 0; i < N; ++i)
					diagD[i] /= sigma;
			}
		} else if (math.max(UBound) < Double.MAX_VALUE && math.min(LBound) > -Double.MAX_VALUE) { // using
			// boundaries
			for (i = 0; i < N; ++i)
				diagD[i] = 0.3 * (UBound[i] - LBound[i]);
			sigma = math.max(diagD);
			if (math.min(diagD) <= 0)
				error("lower boundary must be smaller than upper boundary");
			for (i = 0; i < N; ++i)
				diagD[i] /= sigma;
		} else {
			error("no initial standard deviation or search interval specified");
			sigma = 0.5;
		}

		if (sigma <= 0 || math.min(diagD) <= 0) {
			error("initial standard deviations not specified or non-positive, use setInitialStandarddeviations() or \n" + "    or setInitialSearchRegion().");
			sigma = 1;
		}
		/* save initial standard deviation */
		if (startsigma == null || startsigma.length == 1) {
			startsigma = new double[N];
			for (i = 0; i < N; ++i) {
				startsigma[i] = sigma * diagD[i];
			}
		}
		maxstartsigma = math.max(startsigma);
		minstartsigma = math.min(startsigma);
		axisratio = maxstartsigma / minstartsigma; // axis parallel
		// distribution

		if (xmean == null || sigma == 0)
			error("Bug in init()");

		/* non influencible parameters */
		pc = new double[N];
		ps = new double[N];
		B = new double[N][N];
		C = new double[N][N];
		mC = new Matrix(C);

		xold = new double[N];
		BDz = new double[N];
		bestever = new SolutionGeneric(xmean);
		artmp = new double[N];

		fit.deltaFitHist = new double[5];
		fit.idxDeltaFitHist = -1;
		for (i = 0; i < fit.deltaFitHist.length; ++i)
			fit.deltaFitHist[i] = 1.;

		// code dublicated in reSizeLambda
		fit.fitness = new IntDouble[sp.getLambda()]; // including penalties
		fit.funValues = new IntDouble[sp.getLambda()]; // raw function values
		fit.history = new double[10 + 30 * N / sp.getLambda()];

		arx = new double[sp.getLambda()][N];
		arxrepaired = new double[sp.getLambda()][N];

		for (i = 0; i < sp.getLambda(); ++i) {
			fit.fitness[i] = new IntDouble();
			fit.funValues[i] = new IntDouble();
		}

		// initialization
		for (i = 0; i < N; ++i) {
			pc[i] = 0;
			ps[i] = 0;
			for (int j = 0; j < N; ++j) {
				B[i][j] = 0;
				C[i][j] = 0;
			}
			B[i][i] = 1;
			C[i][i] = diagD[i] * diagD[i];
		}
		maxsqrtdiagC = Math.sqrt(math.max(math.diag(C)));
		minsqrtdiagC = Math.sqrt(math.min(math.diag(C)));
		countCupdatesSinceEigenupdate = 0;
		iniphase = false; // obsolete

		// fileswritten = new String[]{""};
		state = 0; // init done
		timings.start = System.currentTimeMillis();

		/* Some consistency check */
		for (i = 0; i < N; ++i) {
			if (LBound[i] > UBound[i])
				error("lower bound is greater than upper bound");
			if (initXLow != null)
				if (LBound[i] > initXLow[i])
					error("lower bound '" + LBound[i] + "'is greater than lower initial search region" + initXLow[i]);
			if (initXUp != null)
				if (UBound[i] < initXUp[i])
					error("upper bound '" + UBound[i] + "' is smaller than upper initial search region " + initXUp[i]);
		}
		String[] s = stopConditions.getStrings();
		if (!s[0].equals(""))
			warning("termination condition satisfied by initialization: \n  " + s[0]);

	} // init()

	/**
	 * reads and sets properties and options from default input file
	 * "CMAEvolutionStrategy.properties"
	 * 
	 */
	public Properties loadProperties() {
		return loadProperties(propertiesFileName);
	}

	Properties properties = new Properties();

	/**
	 * reads and sets properties and options from default input file fileName
	 * 
	 * @param fileName
	 *            can be empty, where no reading takes place
	 */
	public Properties loadProperties(String fileName) {
		this.propertiesFileName = fileName;
		if (fileName.equals(""))
			return properties;
		try {
			FileInputStream fis = new FileInputStream(fileName);
			properties.load(fis);
			fis.close();
		} catch (IOException e) {
			e.printStackTrace();
			System.out.println("File '" + fileName + "' not found, no options read");
		}
		setProperties(properties);
		return properties;
	}

	/**
	 * set CMA properties and options from Properties class input
	 * 
	 * @param properties
	 *            java.util.Properties key-value hash table
	 * @see #loadProperties()
	 */
	public void setProperties(Properties properties) {
		String s;
		opts.setOptions(properties);
		if (state > 0) // only options can be changed afterwards
			return;
		state = -1; // read new properties
		// [!n] Niko Hansen says new implementation coming up, dont use this
		// one.
		if (properties.containsKey("boundaryLower") && properties.containsKey("boundaryUpper")) {
			setBoundaries(parseDouble(getAllToken(properties.getProperty("boundaryLower"))), parseDouble(getAllToken(properties.getProperty("boundaryUpper"))));
		}

		if (properties.containsKey("initialSearchRegionLower") && properties.containsKey("initialSearchRegionUpper")) {
			setInitialSearchRegion(parseDouble(getAllToken(properties.getProperty("initialSearchRegionLower"))), parseDouble(getAllToken(properties.getProperty("initialSearchRegionUpper"))));
		}
		if ((s = properties.getProperty("initialX")) != null) {
			setInitialX(parseDouble(getAllToken(s)));
		}
		if ((s = properties.getProperty("initialStandardDeviations")) != null) {
			setInitialStandardDeviations(parseDouble(getAllToken(s)));
		}
		if ((s = properties.getProperty("dimension")) != null) { // parseInt
			// does not
			// ignore
			// trailing
			// spaces
			setDimension(Integer.parseInt(getFirstToken(s)));
		}
		if ((s = properties.getProperty("populationSize")) != null) {
			sp.setPopulationSize(Integer.parseInt(getFirstToken(s)));
		}
	}

	/** should become generic with type argument? */
	private String getFirstToken(String s) {
		if (s == null)
			error("input string must be non-null");
		String[] ar = s.split(new String("\\s+"));
		return ar[0];
	}

	/** Returns defaultS if s==null or empty, should become generic */
	private Double getFirstToken(String s, Double def) {
		if (s == null)
			return def;
		String[] ar = s.split("\\s+");
		if (ar[0].equals(""))
			return def;
		return Double.valueOf(ar[0]);
	}

	private String removeComments(String s) {
		int i;
		// remove trailing comments
		i = s.indexOf("#");
		if (i >= 0)
			s = s.substring(0, i);
		i = s.indexOf("!");
		if (i >= 0)
			s = s.substring(0, i);
		i = s.indexOf("%");
		if (i >= 0)
			s = s.substring(0, i);
		i = s.indexOf("//");
		if (i >= 0)
			s = s.substring(0, i);
		return s;
	}

	/** Returns def if s==null or empty, code dublicate, should become generic */
	private Long getFirstToken(String s, Long def) {
		if (s == null)
			return def;
		String[] ar = removeComments(s).split("\\s+");
		if (ar[0].equals(""))
			return def;
		return Long.valueOf(ar[0]);
	}

	private String[] getAllToken(String s) {
		// split w.r.t. white spaces regexp \s+
		return removeComments(s).split("\\s+");
	}

	private double[] parseDouble(String[] ars) {
		double[] ard = new double[ars.length];
		for (int i = 0; i < ars.length; ++i) {
			ard[i] = Double.parseDouble(ars[i]);
		}
		return ard;
	}

	/**
	 * For future use in adaptive population size schemes, what happens with
	 * sp.lambda then? For multistarts this class should be reconstruced,
	 * because initialization of xmean, sigma, C needs to be done as well.
	 * 
	 */
	protected void reSizeLambda() {
		if (!(state == 0 || state == 3))
			error("population size cannot be changed in state " + state);
		int i;

		if (sp.flgLambdaChanged == 0)
			return;

		state = 0.5; // new population size not yet sampled
		fit.fitness = new IntDouble[sp.getLambda()]; // including penalties
		fit.funValues = new IntDouble[sp.getLambda()]; // raw function values

		double[] t = fit.history;
		fit.history = new double[10 + 30 * N / sp.getLambda()];
		for (i = 0; i < Math.min(fit.history.length, t.length); ++i)
			fit.history[i] = t[i];
		for (; i < fit.history.length; ++i)
			fit.history[i] = fit.history[i - 1];

		arx = new double[sp.getLambda()][N];
		arxrepaired = new double[sp.getLambda()][N];

		for (i = 0; i < sp.getLambda(); ++i) {
			fit.fitness[i] = new IntDouble();
			fit.funValues[i] = new IntDouble();
		}

	}

	// private void info(String s) {
	// System.out.println(" CMA-ES info: " + s);
	// }

	private void infoVerbose(String s) {
		System.out.println(" CMA-ES info: " + s);
	}

	private void warning(String s) {
		System.out.println(" CMA-ES warning: " + s);
	}

	private void error(String s) {
		System.out.println(" CMA-ES error: " + s);
		// e.printStackTrace(); // output goes to Sytem.err
		// e.printStackTrace(System.out); // send trace to stdout

		throw new CMAException(" CMA-ES error: " + s);
		// System.exit(-1);
	}

	/** some simple math utilities */
	class MyMath { // implements java.io.Serializable {
		int itest;

		double square(double d) {
			return d * d;
		}

		double prod(double[] ar) {
			double res = 1.0;
			for (int i = 0; i < ar.length; ++i)
				res *= ar[i];
			return res;
		}

		public double median(double ar[]) {
			// need a copy of ar
			double[] ar2 = new double[ar.length];
			for (int i = 0; i < ar.length; ++i)
				ar2[i] = ar[i];
			Arrays.sort(ar2);
			if (ar2.length % 2 == 0)
				return (ar2[ar.length / 2] + ar2[ar.length / 2 - 1]) / 2.;
			else
				return ar2[ar.length / 2];
		}

		/** @return Maximum value of 1-D double array */
		public double max(double ar[]) {
			int i;
			double m;
			m = ar[0];
			for (i = 1; i < ar.length; ++i) {
				if (m < ar[i])
					m = ar[i];
			}
			return m;
		}

		/** @return index of minium value of 1-D double array */
		public int minidx(double ar[]) {
			return minidx(ar, ar.length - 1);
		}

		/**
		 * @return index of minium value of 1-D double array between index 0 and
		 *         maxidx
		 * @param ar
		 *            double[]
		 * @param maxidx
		 *            last index to be considered
		 */
		public int minidx(double[] ar, int maxidx) {
			int i, idx;
			idx = 0;
			for (i = 1; i < maxidx; ++i) {
				if (ar[idx] > ar[i])
					idx = i;
			}
			return idx;
		}

		/**
		 * @return index of minium value of 1-D double array between index 0 and
		 *         maxidx
		 * @param ar
		 *            double[]
		 * @param maxidx
		 *            last index to be considered
		 */
		protected int minidx(IntDouble[] ar, int maxidx) {
			int i, idx;
			idx = 0;
			for (i = 1; i < maxidx; ++i) {
				if (ar[idx].d > ar[i].d)
					idx = i;
			}
			return idx;
		}

		/** @return index of maximum value of 1-D double array */
		public int maxidx(double ar[]) {
			int i, idx;
			idx = 0;
			for (i = 1; i < ar.length; ++i) {
				if (ar[idx] < ar[i])
					idx = i;
			}
			return idx;
		}

		/** @return Minimum value of 1-D double array */
		public double min(double ar[]) {
			int i;
			double m;
			m = ar[0];
			for (i = 1; i < ar.length; ++i) {
				if (m > ar[i])
					m = ar[i];
			}
			return m;
		}

		/**
		 * @return Maximum value of 1-D Object array where the object implements
		 *         Comparator Example: max(Double arx, arx[0])
		 */
		public Object max(Object ar[], Comparator c) {
			int i;
			Object m;
			m = ar[0];
			for (i = 1; i < ar.length; ++i) {
				if (c.compare(m, ar[i]) > 0)
					m = ar[i];
			}
			return m;
		}

		/** @return Maximum value of 1-D IntDouble array */
		public IntDouble max(IntDouble ar[]) {
			int i;
			IntDouble m;
			m = ar[0];
			for (i = 1; i < ar.length; ++i) {
				if (m.compare(m, ar[i]) < 0)
					m = ar[i];
			}
			return m;
		}

		/** @return Minimum value of 1-D IntDouble array */
		public IntDouble min(IntDouble ar[]) {
			int i;
			IntDouble m;
			m = ar[0];
			for (i = 1; i < ar.length; ++i) {
				if (m.compare(m, ar[i]) > 0)
					m = ar[i];
			}
			return m;
		}

		/** @return Minimum value of 1-D Object array defining a Comparator */
		public Object min(Object ar[], Comparator c) {
			int i;
			Object m;
			m = ar[0];
			for (i = 1; i < ar.length; ++i) {
				if (c.compare(m, ar[i]) < 0)
					m = ar[i];
			}
			return m;
		}

		/**
		 * @return Diagonal of an 2-D double array
		 */
		public double[] diag(double ar[][]) {
			int i;
			double[] diag = new double[ar.length];
			for (i = 0; i < ar.length && i < ar[i].length; ++i)
				diag[i] = ar[i][i];
			return diag;
		}

		/**
		 * @return 1-D double array of absolute values of an 1-D double array
		 */
		public double[] abs(double v[]) {
			double res[] = new double[v.length];
			for (int i = 0; i < v.length; ++i)
				res[i] = Math.abs(v[i]);
			return res;
		}
	}

	class Timing {
		Timing() {
		}

		Timing(long birth) {
			this.birth = birth;
		}

		long birth; // time at construction
		long start; // time after init
		long eigendecomposition; // spent time in eigendecomposition
	}

	Timing timings = new Timing(System.currentTimeMillis());

	void eigendecomposition() {
		/* Update B and D, calculate eigendecomposition */
		int i, j;

		// kkk TODO: Introduce flgForceUpdate?
		// 20% is usually better in terms of running *time*
		if (timings.eigendecomposition < 1000 + opts.maximumTimeFractionForEigendecomposition.doubleValue() * (System.currentTimeMillis() - timings.start) && (true && countCupdatesSinceEigenupdate > 0) && (countCupdatesSinceEigenupdate > 1. / sp.getCcov() / N / 5.)) {
			long firsttime = System.currentTimeMillis();
			if (false) { // for debugging
				try {
					PrintWriter out2 = new PrintWriter(new FileWriter("tmpcmaes.dat", countiter > 1));
					out2.println("countiter=" + countiter);
					for (i = 0; i < N; ++i) {
						for (j = 0; j < N; ++j) {
							out2.print(C[i][j] + " ");
						}
						out2.println("");
					}
					out2.println("");
					out2.flush();
				} catch (IOException e) {
					e.printStackTrace();
				}
			} // end for debugging

			// enforce symmetry of C
			for (i = 0; i < N; ++i)
				for (j = 0; j < i; ++j)
					C[j][i] = C[i][j];

			if (1 == 1) {
				EigenvalueDecomposition ed = new Matrix(C).eig();
				B = ed.getV().getArray(); // is not a copy
				try {
					for (i = 0; i < N; ++i)
						diagD[i] = Math.sqrt(ed.getD().get(i, i));
				} catch (RuntimeException e) {
					error("BUG: eigenvalues should never become negative");
				}
				if (math.min(diagD) == 0) // error management is done
					// elsewhere
					axisratio = Double.POSITIVE_INFINITY;
				else
					axisratio = math.max(diagD) / math.min(diagD);
			}
			countCupdatesSinceEigenupdate = 0;
			timings.eigendecomposition += System.currentTimeMillis() - firsttime;
			// ccc test.BD = true;
		} // end Update B and D
	} // eigendecomposition

	/**
	 * Samples the recent search distribution lambda times.
	 * 
	 * @return lambda by dimension array of sampled solutions, where lambda =
	 *         parameters. getPopulationSize()
	 * @see #updateDistribution(double[])
	 * @see #resampleSingle(int)
	 */
	public double[][] samplePopulation() {
		int i, j, iNk;
		double sum;

		if (state < 0)
			init();
		else if (state < 3 && state > 2)
			error("mixing of calls to updateSingle() and samplePopulation() is not possible");
		else
			eigendecomposition(); // latest possibility to generate B and
		// diagD

		reSizeLambda(); // in case lambda was changed externally or multistart

		if (state != 1)
			++countiter;
		state = 1; // can be repeatedly called without problem
		idxRecentOffspring = sp.getLambda() - 1;

		// maximal and minimal standard deviations
		if (opts.lowerStandardDeviations != null && opts.lowerStandardDeviations.length > 0)
			for (i = 0; i < N; ++i) {
				double d = opts.lowerStandardDeviations[Math.min(i, opts.lowerStandardDeviations.length - 1)];
				if (d > sigma * minsqrtdiagC)
					sigma = d / minsqrtdiagC;
			}
		if (opts.upperStandardDeviations != null && opts.upperStandardDeviations.length > 0)
			for (i = 0; i < N; ++i) {
				double d = opts.upperStandardDeviations[Math.min(i, opts.upperStandardDeviations.length - 1)];
				if (d < sigma * maxsqrtdiagC)
					sigma = d / maxsqrtdiagC;
			}

		testAndCorrectNumerics();

		/* sample the distribution */
		for (iNk = 0; iNk < sp.getLambda(); ++iNk) { /*
														 * generate scaled
														 * random vector (D * z)
														 */
			for (i = 0; i < N; ++i) {
				artmp[i] = diagD[i] * rand.nextGaussian();
			}
			/* add mutation (sigma * B * (D*z)) */
			for (i = 0; i < N; ++i) {
				for (j = 0, sum = 0; j < N; ++j)
					sum += B[i][j] * artmp[j];
				arx[iNk][i] = xmean[i] + sigma * sum;
			}
		}

		/* make a (deep enough) copy of arx to allow repair */
		for (i = 0; i < sp.getLambda(); ++i)
			for (j = 0; j < N; ++j)
				arxrepaired[i][j] = Math.min(Math.max(arx[i][j], LBound[j]), UBound[j]);

		return arxrepaired;
	} // end samplePopulation()

	/**
	 * samples a single new search point.
	 * <p>
	 * Here is an example code how to use the method.
	 * 
	 * <pre>
	 * public class MinimalCodeExample {
	 * 	public static void main(String[] args) {
	 * 		CMAEvolutionStrategy cma = new CMAEvolutionStrategy();
	 * 
	 * 		// setup initialization and options
	 * 		int dim = 10;
	 * 		cma.setDimension(dim);
	 * 		cma.setInitialSearchRegion(3, 7); // init x and sigma concurrently 
	 * 		cma.opts.stopTolFun = new Double(1e-8);
	 * 		cma.opts.stopIterations = new Long(10 + 20 * dim * dim);
	 * 
	 * 		// iteration loop
	 * 		while (cma.stopConditions.isFalse()) {
	 * 			// get new solution point
	 * 			double[] x = cma.sampleSingle();
	 * 
	 * 			// compute objective function value on the fly
	 * 			double value = 0.0;
	 * 			for (int i = 0; i &lt; x.length; ++i)
	 * 				value += (i + 1) * x[i] * x[i];
	 * 
	 * 			// update distribution based on computed function value
	 * 			cma.updateSingle(value);
	 * 
	 * 			// output
	 * 			writeToDefaultFilesNew(); // write data to files
	 * 			if (cma.counteval % 100 == 1)
	 * 				System.out.println(cma.printCaption());
	 * 			if (cma.counteval % 9 == 0)
	 * 				System.out.println(cma.print());
	 * 		}
	 * 		System.out.println(&quot;First stop condition in cma: '&quot; + cma.stopConditions.getStrings()[0] + &quot;'&quot;);
	 * 		// double[] cma.getBestSolution().getX(); // is the best evaluated point
	 * 		// double[] cma.getMeanX();               // was not evaluated but may be even better 
	 * 	} // main
	 * } // class
	 * </pre>
	 * 
	 * @return new solution point as double[]
	 * @see #updateSingle(double)
	 * @see #resampleSingle()
	 */
	public double[] sampleSingle() {
		if (state == 1)
			error("call updateSingle() before calling sampleSingle() again");
		if (mode == PARALLEL_MODE)
			error("sampleSingle() and samplePopulation() cannot be mixed up");

		mode = SINGLE_MODE;

		++idxRecentOffspring;

		if (counteval == 0 || idxRecentOffspring == sp.getLambda()) {
			samplePopulation();
			idxRecentOffspring = 0;
		}
		state = 1;

		return arxrepaired[idxRecentOffspring];
	}

	/**
	 * re-generate an (unsatisfactory) solution. After getting a single point
	 * via calling sampleSingle() this point can be sampled anew by
	 * resampleSingle()
	 */
	public double[] resampleSingle() {
		if (state != 1)
			error("call sampleSingle() before calling resampleSingle()");
		return resampleSingle(idxRecentOffspring);
	}

	/**
	 * re-generate the index-th solution. After getting lambda solution points
	 * with samplePopulation() the i-th point, i=0...lambda-1, can be sampled
	 * anew by resampleSingle(i). Remark that after pop=samplePopulation();
	 * x=resampleSingle(i); we have pop[i]==x.
	 * 
	 * @see #samplePopulation()
	 */
	public double[] resampleSingle(int index) {
		int i, j;
		double sum;
		if (state != 1)
			error("call samplePopulation before calling resampleSingle(int index)");

		/* sample the distribution */
		/* generate scaled random vector (D * z) */
		for (i = 0; i < N; ++i) {
			artmp[i] = diagD[i] * rand.nextGaussian();
		}
		/* add mutation (sigma * B * (D*z)) */
		for (i = 0; i < N; ++i) {
			for (j = 0, sum = 0.0; j < N; ++j)
				sum += B[i][j] * artmp[j];
			arx[index][i] = xmean[i] + sigma * sum;
		}
		/* make a (deep enough) copy of arx to allow repair */
		for (j = 0; j < N; ++j)
			arxrepaired[index][j] = Math.min(Math.max(arx[index][j], LBound[j]), UBound[j]);
		return arxrepaired[index];
	}

	public void updateSingle(double functionValue) {
		if (state < 1)
			error("call sampleSingle() first");
		if (state >= 3)
			error("call sampleSingle() first");
		if (idxRecentOffspring >= sp.getLambda()) // better use assert ?
			error("Bug detected");

		mode = SINGLE_MODE;

		++counteval;

		/* transfer/pass input argument */
		fit.funValues[idxRecentOffspring].d = functionValue;
		fit.funValues[idxRecentOffspring].i = idxRecentOffspring;

		recentFunctionValue = functionValue;
		if (idxRecentOffspring == 0) {
			recentMaxFunctionValue = recentFunctionValue;
			recentMinFunctionValue = recentFunctionValue;
		} else {
			recentMaxFunctionValue = Math.max(recentMaxFunctionValue, recentFunctionValue);
			recentMinFunctionValue = Math.min(recentMinFunctionValue, recentFunctionValue);
		}

		state = 2.5; // pretty ugly due to historical reasons

		if (idxRecentOffspring == sp.getLambda() - 1) {
			updateDistribution();
		} else if (functionValue < bestever.getFunctionValue() || counteval == 1) {
			/* assign best solution */
			bestever = new SolutionGeneric(arxrepaired[idxRecentOffspring], functionValue, counteval);
		}
	}

	/**
	 * update of the search distribution after samplePopulation(). Should not be
	 * mixed with calls to sampleSingle() or updateSingle().
	 * 
	 * @see #samplePopulation()
	 */
	public void updateDistribution(double[] functionValues) {
		if (state == 3) {
			error("updateDistribution() was already called");
		}
		if (state > 2)
			warning("in method updateDistribution(): previous call to updateSingle() has no effect");

		if (mode == SINGLE_MODE)
			error("sampleSingle() and updateDistribution cannot be mixed up");

		/* transfer/pass input argument */
		for (int i = 0; i < sp.getLambda(); ++i) {
			fit.funValues[i].d = functionValues[i];
			fit.funValues[i].i = i;
		}

		counteval += sp.getLambda();
		recentFunctionValue = math.min(fit.funValues).d;
		recentMaxFunctionValue = math.max(fit.funValues).d;
		recentMinFunctionValue = math.min(fit.funValues).d;
		updateDistribution();
	}

	private IntDouble[] computePenalties() {
		int i, j, iNk;
		/* penalize repairment, eg. for boundaries */
		// TODO: figure out whether the change of penalty is too large or fast
		// which can disturb selection
		// this depence in particular on the length of fit.medianDeltaFit
		if (true || countiter < fit.deltaFitHist.length || countiter % 1 * (N + 2) == 0) {
			// minimum of differences with distance lambda/2, better the
			// 25%tile?
			// assumes sorted array!!
			int ii = (sp.getLambda()) / 2;
			double medianDeltaFit = Math.abs(fit.funValues[ii].d - fit.funValues[0].d);
			for (i = 1; i + ii < sp.getLambda(); ++i)
				// minimum because of outliers
				medianDeltaFit = Math.min(medianDeltaFit, Math.abs(fit.funValues[ii + i].d - fit.funValues[i].d));
			medianDeltaFit /= sigma * sigma; // should be somehow constant,
			// because dfit depends on sigma
			// (verified on sphere)
			if (medianDeltaFit > 0) {
				// System.out.println("set" + medianDeltaFit + " " +
				// math.median(fit.medianDeltaFit));
				if (fit.idxDeltaFitHist == -1) // first time: write all fields
					for (i = 0; i < fit.deltaFitHist.length; ++i)
						fit.deltaFitHist[i] = medianDeltaFit;
				if (++fit.idxDeltaFitHist == fit.deltaFitHist.length)
					fit.idxDeltaFitHist = 0;
				// save last five values in fit.medianDeltaFit
				fit.deltaFitHist[fit.idxDeltaFitHist] = medianDeltaFit;
			}
		}
		/* calculate fitness by adding function value and repair penalty */
		double penfactor = 1. * 5. * math.median(fit.deltaFitHist);
		for (iNk = 0; iNk < sp.getLambda(); ++iNk) {
			double sqrnorm = 0;
			double prod = Math.pow(math.prod(diagD), 1.0 / (double) N);
			/*
			 * calculate C^-1-norm of Delta x: norm(D^(-1) * B^(-1) * (Delta
			 * x))^2
			 */
			for (i = 0; i < N; ++i) {
				double sum = 0.0;
				for (j = 0, sum = 0.; j < N; ++j)
					sum += B[j][i] * ((arxrepaired[fit.funValues[iNk].i][j] - arx[fit.funValues[iNk].i][j]));
				sqrnorm += math.square(sum / (Math.pow(diagD[i], 0.9) * Math.pow(prod, 0.10))); // regularization
				// to I
			}
			// sqrnorm/N equals approximately 1/sigma^2
			fit.fitness[iNk].d = fit.funValues[iNk].d + penfactor * sqrnorm / (N + 2); // /
			// (sigma
			// *
			// sigma);
			fit.fitness[iNk].i = fit.funValues[iNk].i;
			// System.out.println(math.median(fit.medianDeltaFit) + " " +
			// sqrnorm / (N+2)); // / (sigma * sigma));
		}
		// if (countiter % 10 == 1)
		// System.out.println(math.median(fit.medianDeltaFit) + " " + sqrnorm);
		return fit.fitness;

	}

	private void updateDistribution() {

		int i, j, k, iNk, hsig;
		double sum, tfac;
		double psxps;

		if (state == 3) {
			error("updateDistribution() was already called");
		}

		/* sort function values */
		Arrays.sort(fit.funValues, fit.funValues[0]);

		/*
		 * compute penalized fitness values based on difference between
		 * arxrepaired and arx
		 */
		fit.fitness = computePenalties(); // assignment is superfluous

		/* sort fitness values */
		java.util.Arrays.sort(fit.fitness, fit.fitness[0]);

		/* update fitness history */
		for (i = fit.history.length - 1; i > 0; --i)
			fit.history[i] = fit.history[i - 1];
		fit.history[0] = fit.funValues[0].d;

		/* save/update bestever-value */
		if (countiter == 1 || fit.funValues[0].d < bestever.functionValue) {
			bestever.functionValue = fit.funValues[0].d;
			bestever.evaluation = counteval - sp.getLambda() + fit.funValues[0].i + 1;
			bestever.x = arxrepaired[fit.funValues[0].i];
		}

		/* calculate xmean and rgBDz~N(0,C) */
		for (i = 0; i < N; ++i) {
			xold[i] = xmean[i];
			xmean[i] = 0.;
			for (iNk = 0; iNk < sp.getMu(); ++iNk)
				xmean[i] += sp.getWeights()[iNk] * arx[fit.fitness[iNk].i][i];
			BDz[i] = Math.sqrt(sp.getMueff()) * (xmean[i] - xold[i]) / sigma;
		}

		/* calculate z := D^(-1) * B^(-1) * BDz into artmp */
		for (i = 0; i < N; ++i) {
			for (j = 0, sum = 0.; j < N; ++j)
				sum += B[j][i] * BDz[j];
			artmp[i] = sum / diagD[i];
		}
		/* cumulation for sigma (ps) using B*z */
		for (i = 0; i < N; ++i) {
			for (j = 0, sum = 0.; j < N; ++j)
				sum += B[i][j] * artmp[j];
			ps[i] = (1. - sp.getCs()) * ps[i] + Math.sqrt(sp.getCs() * (2. - sp.getCs())) * sum;
		}

		/* calculate norm(ps)^2 */
		psxps = 0;
		for (i = 0; i < N; ++i)
			psxps += ps[i] * ps[i];

		/* cumulation for covariance matrix (pc) using B*D*z~N(0,C) */
		hsig = 0;
		if (Math.sqrt(psxps) / Math.sqrt(1. - Math.pow(1. - sp.getCs(), 2 * countiter)) / sp.chiN < 1.5 + 1. / (N - 0.5)) {
			hsig = 1;
		}
		for (i = 0; i < N; ++i) {
			pc[i] = (1. - sp.getCc()) * pc[i] + hsig * Math.sqrt(sp.getCc() * (2. - sp.getCc())) * BDz[i];
		}

		/* stop initial phase */
		if (iniphase && countiter > Math.min(1 / sp.getCs(), 1 + N / sp.getMucov()))
			if (psxps / sp.getDamps() / (1. - Math.pow((1. - sp.getCs()), countiter)) < N * 1.05)
				iniphase = false;

		/*
		 * TODO: reconsider this, it is harmful in a dynamic environment remove
		 * momentum in ps, if ps is large and fitness is getting worse
		 */
		if (11 < 3 && psxps / N > 1.5 + 10 * Math.sqrt(2. / N) && fit.history[0] > fit.history[1] && fit.history[0] > fit.history[2]) {
			infoVerbose(countiter + ": remove momentum " + psxps / N + " " + ps[0] + " " + sigma);
			// for (i=1;i<N;++i)
			// System.out.println(ps[i]);
			tfac = Math.sqrt((1 + Math.max(0, Math.log(psxps / N))) * N / psxps);
			for (i = 0; i < N; ++i)
				ps[i] *= tfac;
			psxps *= tfac * tfac;
		}

		/* update of C */
		// sp.ccov = 0;
		if (sp.getCcov() > 0 && iniphase == false) {

			++countCupdatesSinceEigenupdate;

			/* update covariance matrix */
			for (i = 0; i < N; ++i)
				for (j = 0; j <= i; ++j) {
					C[i][j] = (1 - sp.getCcov()) * C[i][j] + sp.getCcov() * (1. / sp.getMucov()) * (pc[i] * pc[j] + (1 - hsig) * sp.getCc() * (2. - sp.getCc()) * C[i][j]);
					for (k = 0; k < sp.getMu(); ++k) { /*
														 * additional rank mu
														 * update
														 */
						C[i][j] += sp.getCcov() * (1 - 1. / sp.getMucov()) * sp.getWeights()[k] * (arx[fit.fitness[k].i][i] - xold[i]) * (arx[fit.fitness[k].i][j] - xold[j]) / sigma / sigma;
					}
				}
			maxsqrtdiagC = Math.sqrt(math.max(math.diag(C)));
			minsqrtdiagC = Math.sqrt(math.min(math.diag(C)));
		} // update of C

		/* update of sigma */
		sigma *= Math.exp(((Math.sqrt(psxps) / sp.chiN) - 1) * sp.getCs() / sp.getDamps());

		state = 3;

	} // updateDistribution()

	/**
	 * runs the optimization for a given number of evaluations on
	 * ObjectiveFunction f to be minimized.
	 * 
	 * @param f
	 *            the ObjectiveFunction to be minimized
	 * @param evaluations
	 *            the int number of evaluations to be executed
	 */
	public void evolveSingle(IObjectiveFunction f, int evaluations) {
		for (int g = 0; stopConditions.test() == 0 && g < evaluations; ++g) {
			updateSingle(f.valueOf(sampleSingle()));
			if (opts.flgWriteFile > 0)
				writeToDefaultFilesNew(); // TODO
			// write
			// depending
			// on time
			// spent on
			// writing
		}

	}

	/**
	 * runs the optimization for a given number of iterations (generations) on
	 * ObjectiveFunctionParallel f to be minimized.
	 * 
	 * <p>
	 * Example code snippet:
	 * 
	 * <pre>
	 * public static void main(String[] args) {
	 * 	ObjectiveFunctionParallel f = new ElliParallel(1e3);
	 * 	CMAEvolutionStrategy cma = new CMAEvolutionStrategy(10);
	 * 	cma.setBoundaries(-10, 10);
	 * 	cma.setInitialSearchRegion(-1, 5);
	 * 	cma.setInitialStandardDeviation(3);
	 * 	cma.opts.stopFunctionValue = new Double(1e-9);
	 * 	while (cma.stopConditions.isFalse()) {
	 * 		if (cma.countiter % 100 == 0)
	 * 			cma.printlnCaption(System.out);
	 * 		cma.iterate(f, 10); // one iteration using function f
	 * 		cma.println(System.out);
	 * 	}
	 * 	System.out.println(cma.stopConditions.getStrings()[0]);
	 * 	System.out.println(&quot;best function value=&quot; + cma.getBestFunctionValue());
	 * }
	 * </pre>
	 * 
	 * @param f
	 *            the ObjectiveFunctionParallel to be minimized
	 * @param iterations
	 *            the int number of iterations to be executed
	 * 
	 */

	public void iterate(IObjectiveFunctionParallel f, int iterations) {
		for (int g = 0; stopConditions.test() == 0 && g < iterations; ++g) {
			updateDistribution(f.valueOf(samplePopulation()));
			if (opts.flgWriteFile > 0)
				writeToDefaultFilesNew(); // TODO
			// write
			// depending
			// on time
			// spent on
			// writing
		}

		/*
		 * makes only sense for large number of iterations double[][] xx = new
		 * double[1][]; xx[0] = xmean; double[] dd = f.f(xx); ++counteval; if
		 * (dd[0] < bestever.functionValue) { // assign best solution
		 * bestever.functionValue = dd[0]; bestever.evaluation = counteval;
		 * bestever.x = xmean; }
		 */
	}

	/**
	 * ratio between length of longest and shortest axis of the distribution
	 * ellipsoid
	 */
	public double getAxisRatio() {
		return axisratio;
	}

	/**
	 * objective function value of best solution found so far.
	 * 
	 * @return objective function value of best solution found so far
	 * @see #getBestSolution()
	 */
	public double getBestFunctionValue() {
		if (state < 0)
			return Double.NaN;
		return bestever.functionValue;
	}

	/**
	 * get best evaluated solution found so far. Remark that the distribution
	 * mean was not evaluated but is expected to have an even better function
	 * value.
	 * <p>
	 * Example: getBestSolution
	 * 
	 * @return best solution (search point) found so far
	 * @see #getBestX()
	 * @see #getMeanX()
	 */
	public ISolution getBestSolution() {
		return bestever;
	}

	/**
	 * best search point of the recent iteration.
	 * 
	 * @return Returns the recentFunctionValue.
	 * @see #getRecentMinFunctionValue()
	 */
	public double[] getBestRecentX() {
		return arxrepaired[math.minidx(fit.funValues, idxRecentOffspring)];
	}

	/**
	 * get best evaluated search point found so far. Remark that the
	 * distribution mean was not evaluated but is expected to have an even
	 * better function value.
	 * 
	 * @return best search point found so far as double[]
	 * @see #getMeanX()
	 */
	public double[] getBestX() {
		if (state < 0)
			return null;
		return bestever.x;
	}

	public int getDimension() {
		return N;
	}

	/**
	 * number of objective function evaluations counted so far
	 */
	public long getCountEval() {
		return counteval;
	}

	/**
	 * number of iterations conducted so far
	 */
	public long getCountIter() {
		return countiter;
	}

	/**
	 * Get mean of the current search distribution. The mean should be regarded
	 * as the best estimator for the global optimimum at the given iteration.
	 * The return value is <em>not</em> a copy. Therefore do not change it,
	 * unless you know what you are doing.
	 * 
	 * @return mean value of the current search distribution
	 * @see #getBestX()
	 */
	public double[] getMeanX() {
		return xmean;
	}

	/**
	 * objective function value of recent solution.
	 * 
	 * @return Returns the recentFunctionValue.
	 */
	public double getRecentFunctionValue() {
		return recentFunctionValue;
	}

	/**
	 * objective function value of the, so far evaluated, worst solution of the
	 * recent iteration.
	 * 
	 * @return Returns the recentMaxFunctionValue.
	 */
	public double getRecentMaxFunctionValue() {
		return recentMaxFunctionValue;
	}

	/**
	 * objective function value of the, so far evaluated, best solution in the
	 * recent iteration (population)
	 * 
	 * @return Returns the recentMinFunctionValue.
	 * @see #bestever
	 */
	public double getRecentMinFunctionValue() {
		return recentMinFunctionValue;
	}

	/** get used random number generator instance */
	public Random getRand() {
		return rand;
	}

	/**
	 * get properties previously read from a property file.
	 * 
	 * @return java.util.Properties key-value hash table
	 * @see #loadProperties()
	 */
	public Properties getProperties() {
		return properties;
	}

	/** @see #setSeed(long) */
	public long getSeed() {
		return seed;
	}

	/**
	 * Set lower and upper boundary in all variables
	 * 
	 * @param xlow
	 * @param xup
	 */
	// [!n] Niko Hansen says new implementation coming up, dont use this one.
	public void setBoundaries(double xlow, double xup) {
		if (_stableMode == true)
			error("boundaries cannot be set (stability issue - you may override this by setting stableMode flag to false, at your own risks!).");

		int len = 1;
		if (N > 0)
			len = N;
		LBound = new double[len];
		UBound = new double[len];
		for (int i = 0; i < len; ++i) {
			LBound[i] = xlow;
			UBound[i] = xup;
		}
	}

	/**/
	/**
	 * sets lower and upper boundaries in all variables.
	 * 
	 * @param xlow
	 *            lower boundary double[], can be 1-D or of length of the number
	 *            of variables (dimension).
	 * @param xup
	 *            see xlow
	 */
	// [!n] Niko Hansen says new implementation coming up, dont use this one.
	public void setBoundaries(double[] xlow, double[] xup) {
		if (xlow == null || xup == null)
			error("boundaries cannot be null");
		if (xlow.length == 1 && xup.length == 1) {
			setBoundaries(xlow[0], xup[0]);
			return;
		}
		if ((N > 0 && (N != xlow.length || N != xup.length)) || (xlow.length != xup.length))
			error("dimensions of boundaries do not match");
		this.LBound = xlow;
		this.UBound = xup;
		N = xlow.length; // changes N only if N was 0
	}

	/**/

	public void setDimension(int n) {
		if (N > 0 && N != n)
			error("dimension cannot be changed anymore or contradicts to initialX");
		N = n;
		sp.init(N);
	}

	public void setInitialSearchRegion(double xlow, double xup) {
		this.initXLow = new double[] { xlow };
		this.initXUp = new double[] { xup };
	}

	/**
	 * sets lower and upper initial value in all variables. Determines the
	 * initial search point, which is uniformely distributed between xlow and
	 * xup, and the initial standard deviations, as 0.3 of the interval, if not
	 * set otherwise.
	 * 
	 * @param xlow
	 *            lower value(s) double[], can be 1-D or of length of the number
	 *            of variables (dimension).
	 * @param xup
	 *            see xlow
	 */
	public void setInitialSearchRegion(double[] xlow, double[] xup) {
		if (xlow == null || xup == null)
			error("input to initial search region cannot be null pointer");
		if (xlow.length == 1 && xup.length == 1) {
			setInitialSearchRegion(xlow[0], xup[0]);
			return;
		}
		if (N > 0 && (N != xlow.length || N != xup.length))
			error("dimensions do not match");
		this.initXLow = xlow;
		this.initXUp = xup;
		N = xlow.length; // changes N only if N was 0
	}

	public void setInitialStandardDeviation(double startsigma) {
		this.startsigma = new double[1];
		this.startsigma[0] = startsigma;
	}

	public void setInitialStandardDeviations(double[] startsigma) {
		// assert startsigma != null;
		if (N > 0 && startsigma.length > 1 && N != startsigma.length)
			error("dimensions do not match");
		if (startsigma.length > 1)
			N = startsigma.length;
		this.startsigma = new double[startsigma.length];
		for (int i = 0; i < startsigma.length; ++i)
			this.startsigma[i] = startsigma[i];
	}

	public void setInitialX(double xmean) {
		this.xmean = new double[1];
		this.xmean[0] = xmean;
	}

	public void setInitialX(double[] xmean) {
		if (N > 0 && xmean.length > 1 && N != xmean.length)
			error("dimensions do not match");
		// copy input xmean
		this.xmean = new double[xmean.length];
		for (int i = 0; i < xmean.length; ++i)
			this.xmean[i] = xmean[i];
		if (xmean.length > 1)
			N = xmean.length; // changes N only if N was 0
	}

	public void setRand(Random rand) {
		this.rand = rand;
	}

	/**
	 * Setter for the seed for the random number generator
	 * java.util.Random(seed). Changing the seed will only take effect, if the
	 * first search point was not yet generated.
	 * 
	 * @param seed
	 *            a long value to initialize java.util.Random(seed)
	 */
	public void setSeed(long seed) {
		if (state > 0)
			warning("setting seed has no effect at this point");
		else {
			this.seed = seed;
			rand.setSeed(seed);
		}
	}

	/**
	 * printing output in a viewable formatting style. The printing
	 * 
	 * <pre>
	 * Iteration,#Fevals: rb Function Value Delta( best ,worst) |idx: Max SD idx: Min SD  | minsigD  sigma Axisratio&quot;);
	 *   164( 8),   1638: 5.5061568003892640e-08 (-4e-08,3e-08) |  0: 3.3e-05  8: 1.5e-05 | 1.4e-05 5.6e-05   2.34 
	 * </pre>
	 * 
	 * shows the output of printCaption() in the first line and in the second
	 * line
	 * <li>164 iteration number
	 * <li>( 8) recently sampled search point in this iteration,
	 * <li>1638: number of function evaluations
	 * <li>5.5061568003892640e-08 objective function value F of the best point
	 * in the recent generation
	 * <li>(-4e-08, difference between the best ever evaluated function value
	 * to F,
	 * <li>3e-08) | difference between the worst function value of the recent
	 * generation to F
	 * <li>0: index of coordinate with largest standard deviation
	 * <li>3.3e-05 respective standard deviation
	 * <li>8: index of coordinate with smallest standard deviation
	 * <li>1.5e-05 | respective standard deviation
	 * <li>index of coordinate with smallest standard deviation: respective
	 * standard deviation
	 * <li>| 1.4e-05 standard deviation in smallest principal axis direction
	 * <li> 5.6e-05 sigma
	 * <li> 2.34 axisratio, ie. quotient between the standard deviations in
	 * largest an smallest principal axis directions, ie. sqare root of the
	 * quotient between largest and smallest eigenvalue of covariance matrix C
	 * 
	 * @see #printCaption()
	 */
	public String print() {
		/*
		 * String.format(Locale.US, " %1$4d(%2$2d): %3$5d ", new Object[]{ new
		 * Long(countiter), new Integer(idxRecentOffspring), new Long(counteval) }) +
		 * String.format(Locale.US, "%1$.16e (%2$+.0e %3$.0e)", new Object[]{
		 * new Double(recentFunctionValue), new Double(getBestFunctionValue() -
		 * recentFunctionValue), new Double(recentMaxFunctionValue -
		 * recentFunctionValue) }) + String.format(Locale.US, "%1$7.2f ", new
		 * Object[]{ new Double(axisratio) }) + String.format(Locale.US,
		 * "%1$2d:%2$8.1e %3$2d:%4$8.1e", new Object[]{ new
		 * Integer(math.minidx(math.diag(C))), new Double(sigma *
		 * Math.sqrt(math.min(math.diag(C)))), new
		 * Integer(math.maxidx(math.diag(C))), new Double(sigma *
		 * Math.sqrt(math.max(math.diag(C)))) })
		 */
		String s;
		if (state < 0)
			s = new String(new PrintfFormat(Locale.US, " %4d").sprintf(countiter) + new PrintfFormat(Locale.US, "(%2d), ").sprintf(0) + new PrintfFormat(Locale.US, "%6.0d: ").sprintf(counteval));
		else
			s = new String(new PrintfFormat(Locale.US, " %4d").sprintf(countiter) + new PrintfFormat(Locale.US, "(%2d), ").sprintf(idxRecentOffspring + 1) + new PrintfFormat(Locale.US, "%6.0d: ").sprintf(counteval) + new PrintfFormat(Locale.US, "%.16e ").sprintf(recentFunctionValue) + new PrintfFormat(Locale.US, "(%+.0e,").sprintf(getBestFunctionValue() - recentFunctionValue) + new PrintfFormat(Locale.US, "%.0e) | ").sprintf(recentMaxFunctionValue - recentFunctionValue)
					+ new PrintfFormat(Locale.US, "%2d:").sprintf(math.maxidx(math.diag(C))) + new PrintfFormat(Locale.US, "%8.1e ").sprintf(sigma * maxsqrtdiagC) + new PrintfFormat(Locale.US, "%2d:").sprintf(math.minidx(math.diag(C))) + new PrintfFormat(Locale.US, "%8.1e ").sprintf(sigma * minsqrtdiagC) + new PrintfFormat(Locale.US, "| %6.1e ").sprintf(sigma * math.min(diagD)) + new PrintfFormat(Locale.US, "%6.1e ").sprintf(sigma) + new PrintfFormat(Locale.US, "%6.2f").sprintf(axisratio)
					+ new PrintfFormat(Locale.US, "   | %4.1f ").sprintf((System.currentTimeMillis() - timings.start) / 1000.) + new PrintfFormat(Locale.US, "%4.1f ").sprintf(timings.eigendecomposition / 1000.));

		// TODO: clean up this writing, via better option?
		if (opts.verbose > 1)
			writeToFile(opts.outputFileNamesPrefix + "disp" + ".dat", s);
		return s;

		/*
		 * return new String( new Long(countiter) + " " + new
		 * Integer(idxRecentOffspring) + " " + new Long(counteval) + " " + new
		 * Double(recentFunctionValue) // + " " + new Double(FunctionValue() -
		 * recentFunctionValue) // + " " + new Double(recentMaxFunctionValue -
		 * recentFunctionValue) + " " + new Double(axisratio) + " " + new
		 * Integer(math.minidx(math.diag(C))) + " " + new Double(sigma *
		 * Math.sqrt(math.min(math.diag(C)))) + " " + new
		 * Integer(math.maxidx(math.diag(C))) + " " + new Double(sigma *
		 * Math.sqrt(math.max(math.diag(C)))) );
		 */
		/*
		 * formatting template String.format(Locale.US, "%1$6.2e %2$+.0e", new
		 * Object[]{ new Double(), new Double() })
		 * 
		 */
		// out.print(math.min(diagD));
		// out.print(" ");
		// new DecimalFormat("0.00E0").format((3.34)) + " " +
		// (cma.fit.fitness[(cma.parameters.getLambda()/2)].d
		// - cma.fit.fitness[0].d) + "," +
		// cma.fit.fitness[cma.parameters.getLambda()-1].d + ") | " +
	}

	/**
	 * see print. Calls out.println(print())
	 * 
	 * @see #print()
	 */
	public void println(PrintStream out) {
		out.println(print());
	}

	/** prints a caption for the printings of method print(). */
	public String printCaption() {
		String s = new String("Iteration,#Fevals: rb Function Value Delta( best ,worst) |idx: Max SD idx: Min SD  | minsigD  sigma Axisratio | time .in eig");
		// 491( 3), 3924: 1.1245467061992267e+00 (-2e-01,4e-01) 9: 7.8e-05 2:
		// 5.0e-02 | 5.9e-03 1.3e-02 660.41
		// TODO: clean up this writing, via better option?
		if (opts.verbose > 1)
			writeToFile(opts.outputFileNamesPrefix + "disp" + ".dat", s);
		return s;
	}

	/** @see #printCaption() */
	public void printlnCaption(PrintStream out) {
		out.println(printCaption());
	}

	/** prints an initial message of CMA-ES optimizer */
	public String printHeader() {
		String s = new String("(" + sp.getMu() + "," + sp.getLambda() + ")-CMA-ES(mu_eff=" + Math.round(100. * sp.getMueff()) / 100. + "), Ver " + versionNumber + ", dimension=" + N);
		return s;

	}

	/** writing output intended to be further processed, eg. for plotting purpose */
	public void write(PrintWriter out) {
		out.println(write());
		out.flush();
	}

	public void writeCaption(PrintWriter out, String commentsymbol) {
		out.println(writeCaption(commentsymbol));
	}

	public String writeCaption(String commentsymbol) {
		return new String(commentsymbol + " counteval functionvalue sigma sigma*sqrt(max(diag(C)))" + "sigma*sqrt(min(diad(C))) axisratio=max(diag(D))/min(diag(D))  min(diag(D))");
	}

	/** writing output intended to be further processed, eg. for plotting purpose */
	public String write() {
		String s = new String("");
		// s = s + (countiter + " ");
		s = s + (counteval + " ");
		s = s + bestever.getFunctionValue();
		// out.print(fit.history[0]);
		s = s + (recentFunctionValue);
		s = s + (" ");
		s = s + (sigma);
		s = s + (" ");
		s = s + (sigma * Math.sqrt(math.max(math.diag(C))));
		s = s + (" ");
		s = s + (sigma * Math.sqrt(math.min(math.diag(C))));
		s = s + (" ");
		s = s + (math.max(diagD) / (math.min(diagD) != 0.0 ? math.min(diagD) : Double.MIN_VALUE));
		s = s + (" ");
		s = s + (math.min(diagD)); // 8th entry
		s = s + (" ");

		double[] tmp = (double[]) diagD.clone();
		java.util.Arrays.sort(tmp);
		for (int i = 0; i < N; ++i) {
			s = s + (tmp[i]);
			s = s + (" ");
		}
		return s;
		// out = new PrintWriter(new FileWriter("xcmaes.dat",
		// countiter>1));
	}

	/**
	 * countiter counteval bestever_function_value recent_function value
	 */
	public String writeSingleFitness() {
		String s = new String();
		s = countiter + " " + counteval + " " + bestever.getFunctionValue() + " " + recentFunctionValue;
		return s;
	}

	public String writeFitness() {
		String s = new String();
		s = countiter + " " + counteval + " " + bestever.getFunctionValue() + " ";
		if (mode == SINGLE_MODE)
			s += recentFunctionValue + " ";
		else {
			s += fit.funValues[0].d + " ";
			s += fit.funValues[sp.getLambda() / 2].d + " ";
			s += fit.funValues[sp.getLambda() - 1].d + " ";
			// for (int i = 0; i < sp.getLambda(); ++i) {
			// s += fit.funValues[i].d + " ";
			// }
		}
		return s;
	}

	public String writeRecentBestX() {
		int idx = 0;
		if (mode == SINGLE_MODE)
			idx = idxRecentOffspring;
		String s = new String();
		s = countiter + " " + counteval + " " + (state == 1 ? Double.NaN : fit.funValues[idx].d) + " ";
		for (int i = 0; i < N; ++i) {
			s += arx[fit.funValues[idx].i][i] + " ";
		}
		return s;
	}

	public String writeMeanX() {
		String s = new String();
		s = countiter + " " + counteval + " 0 ";
		for (int i = 0; i < N; ++i) {
			s += xmean[i] + " ";
		}
		return s;
	}

	/** 4-th to last column are sorted axis lengths axlen */
	public String writeDiagD() {
		String s = new String();
		s = countiter + " " + counteval + " " + sigma + " ";
		double[] tmp = (double[]) diagD.clone();
		java.util.Arrays.sort(tmp);
		for (int i = 0; i < N; ++i) {
			s += tmp[i] + " ";
		}
		return s;
	}

	public String writeSqrtDiagC() {
		String s = new String();
		s = countiter + " " + counteval + " " + sigma + " ";
		for (int i = 0; i < N; ++i) {
			s += sigma * Math.sqrt(C[i][i]) + " ";
		}
		return s;
	}

	/**
	 * correlations and covariances of the search distribution. The first,
	 * '%#'-commented row contains itertation number, evaluation number, and
	 * sigma. In the remaining rows the upper triangular part contains variances
	 * and covariances sigma*sigma*c_ij. The lower part contains correlations
	 * c_ij / sqrt(c_ii * c_jj).
	 */
	public String writeC() {
		int i, j;
		String s = new String();
		s = "%# " + countiter + " " + counteval + " " + sigma + "\n";
		for (i = 0; i < N; ++i) {
			for (j = 0; j < i; ++j)
				// ouput correlation in the lower half
				s += C[i][j] / Math.sqrt(C[i][i] * C[j][j]) + " ";
			for (j = i; j < N; ++j)
				s += sigma * sigma * C[i][j] + " ";
			s += "\n";
		}
		return s;
	}

	private String[] fileswritten = new String[] { "" }; // also

	// (re-)initialized
	// in init()

	/**
	 * writes a string to a file, overwrites first, appends afterwards.
	 * <p>
	 * Example: cma.writeToFile("cmaescorr.dat", cma.writeC());
	 * 
	 * @param filename
	 *            is a String giving the name of the file to be written
	 * @param data
	 *            is a String of text/data to be written
	 */
	public void writeToFile(String filename, String data) {
		boolean appendflag = opts.flgAppendFiles;
		for (int i = 0; !appendflag && i < fileswritten.length; ++i)
			if (filename.equals(fileswritten[i])) {
				appendflag = true;
			}
		PrintWriter out = null;
		try {
			out = new PrintWriter(new FileWriter(filename, appendflag));
			out.println(data);
			out.flush(); // no idea whether this makes sense
			out.close();
		} catch (FileNotFoundException e) {
			warning("Could not find file '" + filename + "'(FileNotFoundException)");
		} catch (IOException e) {
			warning("Could not open/write to file " + filename);
			// e.printStackTrace(); // output goes to Sytem.err
			// e.printStackTrace(System.out); // send trace to stdout
		} finally {
			if (out != null)
				out.close();
		}
		// if first time written
		// append filename to fileswritten
		if (appendflag == false) {
			String s[] = fileswritten;
			fileswritten = new String[fileswritten.length + 1];
			for (int i = 0; i < s.length; ++i)
				fileswritten[i] = s[i];
			fileswritten[fileswritten.length - 1] = new String(filename);
		}
	}

	/**
	 * writes data output. Uses opts.outputFileNamesPrefix to create filenames.
	 * Columns 1-2 are iteration number and function evaluation count, columns
	 * 4- are the data according to the filename.
	 * 
	 * @see #writeToDefaultFilesNew(String fileNamePrefix)
	 */
	public void writeToDefaultFilesNew() {
		writeToDefaultFilesNew(opts.outputFileNamesPrefix);
	}

	/**
	 * Writes data to files <tt>fileNamePrefix</tt>fit.dat, ...xmean.dat
	 * ...xbest.dat, ...std.dat, ...axlen.dat.
	 * 
	 * @see #writeToDefaultFilesNew()
	 * @param fileNamePrefix
	 *            prefix String for filenames created to write data
	 */
	public void writeToDefaultFilesNew(String fileNamePrefix) {
		if (state < 1)
			return;
		writeToFile(fileNamePrefix + "fit.dat", writeFitness());
		writeToFile(fileNamePrefix + "xmean.dat", writeMeanX());
		writeToFile(fileNamePrefix + "xbest.dat", writeRecentBestX());
		writeToFile(fileNamePrefix + "std.dat", writeSqrtDiagC()); // sigma*sqrt(diag(C))
		writeToFile(fileNamePrefix + "axlen.dat", writeDiagD());
	}

	public void writeToDefaultFilesHeaders() {
		writeToDefaultFilesHeaders(opts.outputFileNamesPrefix);
	}

	/**
	 * Writes headers (column annotations) to files <prefix>fit.dat,
	 * ...xmean.dat ...xbest.dat, ...std.dat, ...axlen.dat.
	 * 
	 * @param fileNamePrefix
	 *            String for filenames created to write data
	 */
	public void writeToDefaultFilesHeaders(String fileNamePrefix) {
		writeToFile(fileNamePrefix + "fit.dat", "%# iteration evaluation bestever best median worst\n");
		writeToFile(fileNamePrefix + "xmean.dat", "%# iteration evaluation void mean(1...dimension)\n");
		writeToFile(fileNamePrefix + "xbest.dat", "%# iteration evaluation fitness_of_recent_best x_of_recent_best(1...dimension)\n");
		writeToFile(fileNamePrefix + "std.dat", "%# iteration evaluation sigma sigma*sqrt(diag(C))\n");
		writeToFile(fileNamePrefix + "axlen.dat", "%# iteration evaluation sigma sort(diag(D)) (square roots of eigenvalues of C)\n");
	}

	/**
	 * obsolete, subject to removal
	 */
	public void writeToDefaultFilesOld() {
		int i;
		if (state < 1)
			return;

		// Fitness values

		writeToFile("rescmaes.dat", write());

		// prepare string of xmean to write
		// out = new StringWriter();
		String s = new Long(counteval).toString() + " ";
		for (i = 0; i < N; ++i)
			s = s + new Double(xmean[i]).toString() + " ";

		writeToFile("xcmaes.dat", s);

	}

	/**
	 * very provisorial error handling. Methods of the class
	 * CMAEvolutionStrategy might throw the CMAException, that need not be
	 * catched, because it extends the "unchecked" RuntimeException class
	 */
	public class CMAException extends RuntimeException {
		private static final long serialVersionUID = 1L;

		CMAException(String s) {
			super(s);
		}
	}

	/**
	 * use this method to override the default stable mode (necessary if you
	 * need to set boundaries)
	 * 
	 * @param __mode
	 */
	public void setStableMode(boolean __mode) { // [!n] override boundaries
		// limitation
		this._stableMode = __mode;
		System.out.println("[Warning] overriding stable mode is at your own risks. Boundaries implementation is cma-es may not be optimal.");
	}

}

class IntDouble implements Comparator {
	int i; // unique integer value, useful after sorting
	double d; // double value

	public IntDouble(double d, int i) {
		this.d = d;
		this.i = i;
	}

	public IntDouble(double d) {
		this.d = d;
	}

	public IntDouble() {
	}

	public int compare(Object in1, Object in2) {
		IntDouble o1 = (IntDouble) in1;
		IntDouble o2 = (IntDouble) in2;
		if (o1.d < o2.d)
			return -1;
		if (o1.d > o2.d)
			return 1;
		if (o1.i < o2.i)
			return -1;
		if (o1.i > o2.i)
			return 1;
		return 0;
	}

	public boolean equals(IntDouble o1, IntDouble o2) {
		if (o1.compare(o1, o2) == 0) // && o1.hashCode() == o2.hashCode()
			return true;
		return false;
	}
} // IntDouble
